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A new set of equations for relativistic viscous hydrodynamics that captures both 
On ' weak-coupling and strong-coupling physics to second order in gradients has been 

■ developed recently. We apply this framework to bulk physics at RHIC, both for 
CN ■ standard (Glauber-type) as well as for Color-Glass-Condensate initial conditions 

^ , and show that the results do not depend strongly on the values for the second- 

\ order transport coefficients. Results for multiplicity, radial flow and elliptic flow 

are presented and we quote the ratio of viscosity over entropy density for which 

■ our hydrodynamic model is consistent with experimental data. For Color-Glass- 
Condensate initial conditions, early thermalization does not seem to be required in 

order for hydrodynamics to describe charged hadron elliptic flow. 
-4— • , 

^ ■ I. INTRODUCTION 
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The experimental program at the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven 
has generated a wealth of data B a a i on QCD matter at the highest energy densities 
yjn \ obtained in the laboratory. Remarkably, ideal hydrodynamics seems to offer a sensible 
description of the experimental data for bulk properties (multiplicity, radial and elliptic 
flow) of low pt particles for heavy-ion collisions at RHIC [1, a 0, a 0] ■ 

Upon closer inspection, however, not all of this success can be attributed to modeling 



O 



O ' the system as an ideal fluid. For instance, the energy density distribution which is used 



O 



as an initial condition for the hydrodynamic equations is customarily chosen such that the 
output from the hydrodynamic model matches the experimental data for the multiplicity. 
Furthermore, the time where the hydrodynamic model is initialized as well as the tempera- 
^ ' ture (or energy density) at which the hydrodynamic evolution is stopped are typically chosen 

■ such that the model output matches the experimental data for the radial flow. After these 
parameters have been fixed, only the good description of experimental data for the elliptic 
flow coefficient can be considered a success for ideal hydrodynamics (in the sense that it is 
parameter- free) . 

In order to make progress and learn more about the properties of matter created at RHIC, 
the task is now to both test and improve this ideal hydrodynamic model. The obvious 
framework for this task is dissipative hydrodynamics, since it contains ideal hydrodynamics 
as the special case when all dissipative transport coefficients (such as shear and bulk viscosity 
and heat conductivity) are sent to zero. If the value of the transport coefficients were known 
(e.g. by some first principle calculation), then one could use dissipative hydrodynamics 
to constrain e.g. the initial energy density distribution, which is chosen conveniently in 
the ideal hydrodynamic models. Or otherwise, choosing again physically acceptable initial 
conditions, one is able to constrain the allowed ranges of the transport coefficients. Despite 
recent progress in first principle calculations 11, H, 3, 14, H, 3, 17, 1^, 1^, the values 
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of the hydro dynamic transport coefficients for QCD in the relevant energy range are poorly 
constrained to date, so the second option is currently the only viable possibility. 

For RHIC, the ffist step in this direction was carried out by Teaney [20,], who provided 
estimates for the sign and size of corrections due to shear viscosity. This famous calculation, 
however, did not provide a description of experimental data for non-zero viscosity, because 
it was not dynamic and the initial conditions could not be altered. Only very recently, the 
first hydrodynamic calculations with shear viscosity describing particle spectra for central 
and non-central collisions at RHIC have became available 
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Several other groups haye produced numerical codes capable of performing similar match- 



ing to data [2J, l25|, M, l27|, l28|, l29|, l30|, l31 



However, the precise formulation of the viscous hydrodynamic equations themselves has 
long been debated. To appreciate the complication, one first has to understand a hydrody- 
namic formulation for RHIC physics necessarily has to be fully relativistic, and that the rel- 
ativistic generalization of the Navier-Stokes equations are acausal since they contain modes 
that transport information at superluminal speeds. These are high wavenumber modes and 
therefore in principle outside the range of validity of hydrodynamics, but in practice one has 
to find a way to deal with them in viscous hydrodynamic simulations. A possible solution to 
this problem is known as Miiller-Israel-Stewart theory, where for each transport coefficient a 
corresponding relaxation time is introduced which controls the speed of signal propagation 
for the high wavenumber modes For low momentum modes (up to ffist order 

in gradients), the Miiller-Israel-Stewart theory is identical to the Navier-Stokes equations, 
but differs for higher order gradients. Unfortunately, this implied that the resulting equa- 
tions retained a certain degree of arbitrariness, as it was not clear which additional terms of 
second or higher order in gradients either within the Miiller-Israel-Stewart or other frame- 
works 



[see e.g. 



37, 



1^1) were allowed. For the case of non- vanishing shear viscosity 
only, it was shown recently [40] that the most general form implies five independent terms of 
second order in gradients. This form is general enough to describe the hydrodynamic prop- 
erties of (conformal) plasmas both for weakly coupled systems describable by the Boltzmann 
equation as well as infinitely strongly coupled plasmas, which are accessible via Maldacena's 
conjecture (4l| . 

The aim of this work is to now apply this new set of equations for relativistic shear viscous 
hydrodynamics to the problem of heavy-ion collisions at RHIC. In section [Tll we review the 
setup of conformal relativistic viscous hydrodynamics and our numerics for the simulation 
of heavy-ion collisions. In section UTTt details about the two main models of initial conditions 
for hydrodynamics are given. Section HVl contains our results for the multiplicity, radial flow 
and elliptic flow in Au+Au collisions at top RHIC energies, as well as a note on the notion 
of "early thermalization" . We conclude in section jVl 



II. SETUP 

The energy-momentum tensor for relativistic hydrodynamics in the presence of shear 
viscosity can be written as 

T^'^ = eu^'u^ - pA^"" + H'^^ ( 1 ) 

where e and p are the energy density and pressure, related by an equation of state p = p(e). 

is the fluid four- velocity which fulfills g^^u^^u'^ = 1, where the signature of the metric 
is g^i, = (+, — , — , — ). The projector A^'^ = g^" — u^u^ is orthogonal to the fluid velocity 
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u^A'^" = 0. U^'^ is the viscous shear tensor which is symmetric, traceless (11^ = 0) and 
orthogonal to the fluid velocity. Hydrodynamics describes the evolution of the energy density 
and fluid velocity. The evolution equations are simply given by the conservation of the energy 
momentum tensor D^T^"" = 0, where is the (geometric) covariant derivative. Projection 
of Ua and on D^T^"^ = gives 

De = -{e + p)VX + ^n'^'V(,M^) , (2) 

where D = u"'Da and = A^"Dq, can be thought of as comoving time and space deriva- 
tives, respectively. Note that Z)^ = Ufj,D + V^. The brackets ( ) denote the combination 

Ai,B^) = (a^A^ + A^A^ - ^A"^A,.) A^B,, (3) 

which is a projector that is symmetric, traceless, and orthogonal to the fluid velocity. For 
later convenience, we also introduce symmetric and anti-symmetric brackets 

Ai^^B,) = ^{A^B, + A,B^) 

A[^B,] = ^ (A^B, - A,B^) . (4) 

The equations ([2]) can be considered four equations for the four independent components 
of e,u^. A theory of viscous hydrodynamics still has to specify the evolution or defining 
equations for the five independent components of the shear tensor 11^'^. To first order in 
gradients, these are given by the relativistic Navier-Stokes equations 

W = 7]V^''u^'\ (5) 

where rj is the shear viscosity coefficient. As mentioned in the introduction, this theory 
suffers from acausal signal propagation and associated numerical instabilities. To second 
order in gradients, the evolution equations are given by [13] (see also (i^ l) 

= T^V^^u'^^ - rn [a^A^DH"^ + -n^^(V„M° 



n<^n">^ + ^Y[<\uj'>^ - h^<\uj->^ , (6) 



+ 2 



2r/2 " 2r7 2 

where = — Vj^m,^] is the fluid vorticity and R°'^'^^, R^^ are the Riemann and Ricci tensors, 
respectively. The coefficients rn, /t, Ai, A2, A3 are the five new coefficients controlling the size 
of the allowed terms of second order in gradients. Having an application to the problem of 
heavy-ion collisions in mind, the above set of equations can be simplified: for all practical 
purposes spacetime can be considered flat, such that both the Riemann and Ricci tensors 
vanish identically. Thus, only the four coefficients rn, Ai, A2, A3 enter the problem. 
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A. A note on bulk viscosity and conformality 

Besides shear viscosity, QCD also has no n- vanishing bulk viscosity ( which can be related 
to the QCD trace anomaly ^] 

C ~ = e - 3p. (7) 

QCD lattice simulations seem to indicate that the ratio bulk viscosity over entropy density 
s, (/s, is small compared to rj/s except for a small region around the QCD deconfinement 



transition temperature, where it is sharply peaked |45l . |46| . |47] . If we are interested in 
describing effects from shear viscosity only, we are led to consider ^ = 0, or conformal fluids. 
This has been the main guiding principle in Ref. jiO] and as a consequence Eq. ([6]) obeys 
conformal invariance, unlike most other second-order theories^. 



B. First steps: 0+1 dimensions 

In order to get a crude estimate of the effect of viscous corrections, let us consider the 
arguably simplest model of a heavy-ion collision: a system expanding in a boost-invariant 
fashion along the longitudinal direction and having uniform energy density in the transverse 
plane. Introducing the Milne variables proper time r = yW^^z^ and space-time rapidity 
^ = arctanh(z/t), boost invariance simply translates to requiring all hydrodynamic vari- 
ables (e, M^, n'^'^) to be independent of rapidity, and tensor components m^,!!'^^ to vanish. 
Assuming uniformity in the transverse plane furthermore requires independence from the 
transverse coordinates = {x,y). Even though this means that all the velocity compo- 
nents except u'^ are zero, the system is nevertheless non-trivial in the sense that the sum 
over velocity gradients does not vanish, V^u'^ = ^, sometimes referred to as "Bjorken flow". 

In a way one has modeled an expanding system in static space-time by a system at rest in 
an expanding spacetime. This has been achieved by transforming to the Milne coordinates 
r, ^, where the metric is g^i, = dia.g{grT, dxx, Qyy, 9^^) = (1, — 1, — — t^)- Note that even 
though the spacetime in these coordinates is expanding, it is nevertheless fiat (e.g. has 
vanishing Riemann tensor). 

In this 0+1 dimensional toy model, the viscous hydrodynamic equations become excep- 
tionally simple (4o| . 

e + p Hi 
Ore = -— ^ + 

T T 



m = -'-? + ^ - - TT^ i'^lf ■ (8) 



Tn 3rnr 3r ^ 2rn?7 ^ 

The Navier-Stokes equations are recovered formally in the limit where all second-order co- 
efficients vanish (e.g. rn, Ai -^0); then, one simply has 



The equations ([8]) can be solved numerically along the lines of |38l. |48[|. At very early times. 



where n| > (e + p), the Navier-Stokes equations indicate an increase in energy density and 



^ Note that Muronga derived a version of Eq. ^ in Ref. |36|] that turns out to obey conformal symmetry. 
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a negative effective longitudinal pressure p — n|. Since gradients V ^u^ = 1/r are strongest 
at early times, this suggests that one is applying the Navier-Stokes equations outside their 
regime of validity. Theories including second order gradients may be better behaved at early 
times, but eventually also have to break down when gradients become too strong. Here we 
want to study the effects of the second order coefficients on the value of the shear tensor at 
late times, where a hydrodynamic approach should be valid. 

To this end, let us study the deviation of the shear tensor from its first order value, 
5n = — 1^. At late times, Eq. ([8]) implies e ~ r""^/^, so ?7 ~ r^^. Thus, if 511 is small 
compared to the first order value, from Eq. ([8]) we find 

For a strongly coupled A/" = 4 plasma [13, 113, 42, 43|, one has^ 



^ = -, rn = ^^, Ai = ^, (11) 
s 47r' " 27rT ' ^ 27rT' ^ ' 

and thus n| is larger than its first order value by a factor of 1+ g"!^^ . For RHIC, Tr ^ 1 is 

a reasonable estimate, so one finds that the second order corrections to n| increase its value 
by a few percent over the first order result. 

As an example on the importance of obeying conformal invariance, imagine dropping the 
term involving Vq,m° in the first line of Eq. ([6]). Redoing the above calculation one finds 

which indicates a nearly ten-fold increase of the size of 511 for the non-conformal theory. 



For a weakly coupled plasma well described by the Boltzmann equation [40| , where one has 
Tn = (Ai is unknown but generally set to zero in Miiller-Israel-Stewart theory ), the 
effect may be less pronounced, but still one qualitatively expects second-order effects to be 
anomalously large if conformal invariance is broken in an "ad-hoc" manner. 

Clearly, the above estimates are not meant to be quantitative. Indeed, even the sign of the 
correction may change when allowing more complicated (e.g. three-dimensional) dynamics. 
However, the lesson to be learned from this exercise is that second-order gradients can 
and indeed do modify the shear tensor from its first order (Navier-Stokes) value. This 
is physically acceptable, as long as the second-order corrections are small compared to 
the first order ones (otherwise the system is probably too far from equilibrium for even a 
hydrodynamic description correct to second order in gradients to be valid). A practical 
means for testing this is calculating physical observables for different values of the second- 
order coefficients and making sure that the results do not strongly depend on the choice for 
these specific values. 



C. Including radial flow: lessons from dimensions 

Some more insight on the effect of viscous corrections may be gained by improving the 
model of the previous subsection to allow for radially symmetric dynamics in the transverse 



^ For completeness, we also mention the results n — -pj^,^2 = — , A3 = from 
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plane (but still assuming boost invariance). This is most easily implemented by changing 
to polar coordinates {x,y) — >• (r, 0) with r = v^i?~+^ and = arctan(?//x). In this case, 
the only non- vanishing velocity components are u'^ and , and hence the vorticity u^^" 
vanishes identically. Although non-trivial, the radially symmetric flow case is still a major 
simplification over the general form Eq. ([6]), since again the terms involving k, A2, A3 drop 
out. 

Such a formulation allows both for important code tests [49;] as well as realistic simula- 



tions of central heavy-ion collisions [2]J (note that truncated versions of Eq. ([6]) were used 
in these works). The advantage of this formulation is that since the equations are compar- 
atively simple, it is rather straightforward to implement them numerically and they are not 
very time consuming to solve since only one dimensional (radial) dynamics is involved. The 
shortcoming of simulations with radially symmetric fiow profiles ("radial fiow") is that by 
construction they cannot be matched to experimental data on the impact-parameter depen- 
dence of multiplicity. Thus, the considerable freedom in the initial/final conditions inherent 
to all hydrodynamic approaches cannot be eliminated in this case. 

For this reason, we will choose not to discuss the case of radial fiow here in more detail, 
but rather will comment on it later as a special case of the more general situation. 



D. Elliptic flow: 2-i-l dimensional dynamics 

Retaining the assumption of boost invariance, but allowing for general dynamics in the 
transverse plane, it is useful to keep Cartesian coordinates in the transverse plane, and 
thus u'^,u^,uy are the non-vanishing fiuid velocities. The main reason is that e.g. in polar 
coordinates the equations for the three independent components of U'^'^ would involve some 
extra non- vanishing Christoffel symbols (other than F^^ = t,T^^ = l/r). 

Fortunately, the case of two dimensions is special insofar as the only nontrivial component 
of the vorticity tensor, namely u^^, fulfills the equation [2^ 



xy 



^ „ Dp Du' 

e + p 



om. (13) 



which can be derived by forming the combination V^Du"^ — V^Du^. The expression 0{ll^) 
denotes that the r.h.s. of Eq. f|T3l) is of third order in gradients, and thus should be suppressed 
in the domain of applicability of hydrodynamics. For heavy-ion collisions, typically V ^u^ > 
^, so that for an equation of state with a speed of sound squared = -^M- ~ 4, Eq. ^ 



de 3' 

translates to < unless Dlnu^ > (1 — c^)V^m'^. In particular, this implies that in 

general = is a stable fix-point of the above equation and hence we expect u)^^ to 
remain small throughout the entire viscous hydrodynamic evolution if it is small initially. 

Generically, one uses u^'^ = as an initial condition for hydrodynamics [50|, which implies 
^xy _ g initially. Therefore, to very good approximation we can neglect the terms involving 
vorticity in Eq. ([H]), such that again only the second-order coefficients rn, Ai have to be 
specified. 

The equations to be solved for 2+1 dimensional relativistic viscous hydrodynamics are 
then (in components) 



(e + p)Du' = cl (g'^dje - M*M"9„e) - A^D^H' 
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Dt 



drU" + (9jM* H 

r 



A^'S.uJ' + A"^9,u^' + A^^iM^^ + A'^^iU^ - -A^^/y^u'^ 

3 



-r^A««V, 



(14) 



Here and in the following Latin indices collectively denote the transverse coordinates x, y and 
the relation u^Yi^^ = has been used to derive the above equations (similarly, u^'V(^u^) = 
can be used to obtain the other non-trivial components needed). Note that this particular 
form of Eq. 0141) has not been simplified further since it roughly corresponds to the equations 
implemented for the numerics of and is meant to facilitate understanding of the code 



51[. A simple algorithm to solve Eq. f[T4l) has been outlined in |49| and will be reviewed in 



the next subsection for completeness. 



E. A numerical algorithm to solve relativistic viscous hydrodynamics 



The first step of the algorithm consists of choosing the independent degrees of freedom. 
For boost-invariant 2+1 dimensional dynamics, a sensible choice for this set is e.g. e, u^, u^, 
]jxx^ ]jxy^ ]jyy The pressure is then obtained via the equation of state p{e), and the only 
other non- vanishing velocity as = i/l + + m^. Similarly, the other nonzero components 



of U^'^ are calculated using the equations U'^ = 0, UfJJ^'^ = 0. 

Given the value of the set of independent components at some time r = tq, the aim is 
then to construct an algorithm from Eq. f|T^ such that the new values of the set can be 
calculated as time progresses. Note that in Eq. f|T4|) . time derivatives of the independent 
component set enter only linearly. Therefore, Eq. f|T4l) may be written as a matrix equation 
for the derivatives of the independent component set. 



/ ctoo (^01 ■ 
aio ctii . 



V '^50 0,51 ■ ■ ■ 0,55 / 



«05 \ / \ 
ai5 _ drU^ 



hi 



(15) 



Denoting the above matrix and vector as A and b, respectively, a straightforward way to 
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obtain the time derivatives is via numerical matrix inversion, 



Choosing a naive discretization of derivatives 



A 



-1 



(16) 



drfir) 



fjr + St) - /(r) 
6t 



f{x + a) — f{x — a) 



2a 



(17) 



which is first order accurate in the temporal grid spacing 6t and second order accurate in 
the spatial grid spacing a, one can then directly calculate the new values of the independent 
component set from Eq. (fT6l) . 

Note that for ideal hydrodynamics, the algorithm Eq. (1161) would fail for this naive dis- 



cretization 52 



The reason is that ideal hydrodynamics is inherently unstable to high 
wavenumber fluctuations (which can be thought of as the basis for turbulence). For ideal 
hydrodynamics, one thus has to use a discretization which amounts to the introduction of 
numerical viscosity to dampen these fluctuations. Luckily, viscous hydrodynamics does not 
suffer from this problem because it has real, physical viscosity inbuilt. It is because of this 
reason that the naive discretization can be used in the algorithm Eq. f|T6l) without encoun- 
tering the same problems as in ideal hydrodynamics, as long as a finite value for the viscosity 
7] is used^. While applicable to sufficiently smooth initial conditions, the above algorithm 
is too simple to treat strong gradients such as the propagation of shocks, and should be 
replaced by a more involved scheme in such cases. 

Since matrix inversions are computationally intensive, one can speed up the numerics by 
expressing drll'^" in terms of drU"^, Ore. Inserting these in the equations for Du^ and De, one 
only needs to invert a 3 x 3 matrix to obtain the new values of the energy density and fluid 
velocities. This approach has been used in |2l|, H . 



F. Initial conditions and equation of state 

As outlined in the introduction, any hydrodynamic description of RHIC physics relies on 
given initial energy density distributions. Two main classes of models for boost-invariant 
setups exist: the Glauber models and the Color-Glass-Condensate (CGC) models. 

As will be shown in the following, both model classes can give a reasonable description 
of the experimentally found multiplicity distribution, but they differ by their initial spatial 
eccentricity. A detailed discussion of the initial conditions will be given in subsequent 
sections. 

Besides an initial condition for the energy density, one also needs to specify an initial 
condition for the independent components of the fluid velocities and the shear tensor. For 
the fluid velocities we will follow the standard assumption that these vanish initially (50| . 
Finally, when using the set of equations ([T^ . one also has to provide initial values for the 



^ In practice we have used ^ > 10 ^. Typically, between ^ = 10 ^ and ^ = 10 ^ there are no significant 
changes in the hydrodynamic results and we refer to ^ = 10^'* as "ideal hydrodynamics" . 
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independent components of U^'^. Extreme choices are n^*^ = and a shear tensor so large 
that a diagonal component of the energy-momentum tensor vanishes in the local rest frame 
(e.g. n| = p, or zero longitudinal effective pressure), with the physical result expected 
somewhere in between (see e.g. the discussion in |i53i] ) . 

Once the initial conditions for the independent hydrodynamic variables have been speci- 
fied, one needs the equation of state to solve the hydrodynamic equations ( |T4l) . Aiming for a 
description of deconfined nuclear matter at zero chemical potential, a semi-realistic equation 
of state has to incorporate evidence from lattice QCD calculations [s^] that the transition 
from hadronic to deconfined quark matter is probably an analytic crossover, not a first or 
second order phase transition as often used in ideal hydrodynamic simulations. On the other 
hand, continuum extrapolations for the value of the energy density and pressure for physical 
quark masses are still not accessible with high precision using current lattice methods. For 
this reason, we will employ the equation of state by Laine and Schroder which is derived 
from a hadron resonance gas at low temperatures, a high-order weak-coupling perturbative 
QCD calculation at high temperatures, and an analytic crossover regime interpolating be- 
tween the high and low temperature regime, respectively. For hydrodynamics, an important 
quantity is the speed of sound squared extracted from the equation of state, = For 
completeness, we reproduce a plot of this quantity in Fig. [TJ 



G. Freeze-out 

At some stage in the evolution of the matter produced in a heavy-ion collision, the system 
will become too dilute for a hydrodynamic description to be applicable. This "freeze-out" 
process is most probably happening gradually, but difficult to model realistically. A widely- 
used approximation is therefore to assume instantaneous freeze-out whenever a certain fluid 



cell cools below a certain predefined temperature or energy density (see [30, 157[ for differ- 
ent approaches). The standard prescription for this freeze-out process is the Cooper- Frye 
formula [1^, which allows conversion of the hydrodynamic variables (energy density, fiuid 
velocity,...) into particle distributions. 

Specifically, in the case of isothermal freeze-out at a temperature Tf, the conversion 
from hydrodynamic to particle degrees of freedom will have to take place on a three- 
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dimensional freeze-out hypersurface S, which can be characterized by its normal four-vector, 
and parametrized by three space-time variables jHs], 5^. The spectrum for a single particle 
on mass shell with four momentum p^^ = {E, p) and degeneracy d is then given by 



d 



d^p (27r) 



(18) 



where rfS^ is the normal vector on the hypersurface S and / is the off-equilibrium distribu- 
tion function. 

Originally, the Cooper-Frye prescription was derived for systems in thermal equilibrium, 
where / is built out of a Bose or Fermi distribution function fo{x) = (exp[(x) ±1]~^), 
depending on the statistics of the particle under consideration. In order to generalize it to 
systems out of equilibrium, one customarily relies on the ansatz used in the derivation of 
viscous hydrodynamics from kinetic theory [60| . 



fix^,pn = fo{^) + fo{ 



T 



T 



Mm 



2T^{e + p)' 



(19) 



For simplicity, in the following we approximate fo{x) ~ exp(— x), so similarly 



/(x^p^)=exp(-pX/^) 



1 + 



2T2(e + p) 



(20) 



The effect of this approximation will be commented on in the following sections. 

In practice, for boost-invariant 2+1 dimensional hydrodynamics, the freeze-out hyper- 
surface S'" = (S*, S"^, E^, E^) = (t,x,y,z) can be parametrized either by r, ^ and the polar 
angle 0, or by x,y,^: 



t = T cosh ^ 
X = x(r, 0) 

y = y{'r,4>) ' 

z = T sinh ^ 
The normal vector on E'^ is calculated by 

aE" aE^ dTP 



t = t{x, y) cosh^ 

X = X 

y = y 

z = t{x, y) sinh^ 



(21) 



dE/,(r,0,O = e 



Or 



-drdcpdC, 



and similarly for the other parametrization [61[. 

For a realistic equation of state, at early times the freeze-out hypersurface will contain 
the same transverse coordinate values (x, y) for different times r (see Fig. [2]). Therefore, the 
parametrization in terms of {x,y,^) cannot be used for early times. On the other hand, the 
parametrization in terms of (r, (p, ^) contains derivatives of {x, y) with respect to r, which 
become very large at late times (see Fig. [2]). Numerically, it is therefore not advisable to 
use this parametrization at late times. As a consequence, we use the one parametrization 
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FIG. 2: (Color online) Space-time cut through the three-dimensional hypersurface for a central 
collision within the Glauber model. Simulation parameters used were a = 1 GeV~^, tq = 1 fm/c, 
Tj = 0.36 GeV, Tj = 0.15 GeV, rn = 6^ and Ai = (see next sections for definitions). As can be 
seen from the figure, inclusion of viscosity only slightly changes the form of the surface. 



at early times but switch to the other parametrization at late times, such that the integral 
in Eq. (|T8|) is always defined and numerically well-behaved'^. 

In order to evaluate the integral f|T8|) . it is useful to express also in Milne coordinates, 

p^' = {p\p\py,p^) = (mT cosh(y-0,P^P^ — sinh(F-0), (22) 



where = ym^ + pI + = y E'^ ~ p1- Here and in the following Y = arctanh(p^/i?) 
is the rapidity, and m is the rest mass of the particle under consideration. Then the ^ 
integration can be carried out analytically using 

1 

- / rfecosh"(y-0exp[-;2cosh(F-0] = (-l)"9^i^o(^) = K{n,z), (23) 

where Ko{z) is a modified Bessel function. One finds 

d^N 2d r 

= T^,j dTdct>eM{p^u^+pyuy)/T] X 



mT (^l^ll - {TiK{l, rriTU^/T) + T^^, mru^ /T) + T,K{?,, mTv7 /T)) 

- ( p"It - P^^It I (7^1^(0, rriTU^/T) + T^Kil, rriTU^ /T) + T^K{2, uitU^ /T)) 



Ti = 1 + 



mm + pIw + plwy + 2p^pyWy 



2T\t + p) 



^ It may be possible that other parametrizations may turn out to be more convenient. For instance, it 
is conceivable that performing a triangulation of the three-dimensional hypersurface and replacing the 
integral in (fTS]) by a sum over triangles could turn out to be numerically superior to our method. 
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T2 



-2m' 



2T2(e + p) 



2T2(e + p) 



(24) 



for the {T,(f),C,) parametrization, and a similar result for the other parametrization of the 
hypersurface. The remaining integrals for the particle spectrum have to be carried out 
numerically unless one is considering the case of a central collision [21 , lii where the integral 
has an additional symmetry in 0. 

For the simulation of a heavy-ion collision, one then also needs to take into account the 
feed-down process of particle resonances that decay into lighter, stable particles [gI, 63 . 
Therefore, we calculate the spectra for particle resonances with masses up to ~ 2 GeV and 



then use available routines from the AZHYDRO package [6J] to determine the spectra of 
stable particles including these feed-down contributions. Ultimately, one would be interested 
in describing the last stage of the evolution by coupling the hydrodynamics to a hadronic 



cascade code [65|, |66|, 167|, 168|. We leave this for future work. 

The particle spectra E^^^^ including feed-down contributions can then be used to cal- 
culate experimental observables at central rapidity Y = , such as radial and elliptic flow 
coefficients, fo,f2, respectively. These are defined as 



where 



27r 



E- 



Vo{pt, h) [1 + 2^72 (Pr, h) cos(2(/)p) + 



(25) 



arctan(p^/p^') and px = Jpl + pl- Furthermore, the total multiplicity per unit 



(26) 



rapidity ^ and the mean transverse momentum < pr > are then given by 

dN_^ f, , ^ jdpTplvoipxM 

— = 2tt dpT Pt Vo{pt, 6) , < Pt >= yI 1 m • 

dy J j dpT Pt Vo{pt,o) 

The Pt integrated elliptic fiow coefficient is defined as 



<(6) 



/ dpT PtV2{pt, b)vo{pT, b) 



JdpT Pt Vo{PT,b) 
and the minimum bias elliptic fiow coefficient as 0] 

Jdbb V2{pT,b) vo{pT,b) 



vt{PT) 



Jdbb Vo{pT,b) 



(27) 



(28) 



H. Code tests 

It is imperative to subject the numerical implementation of the relativistic viscous hy- 
drodynamic model to several tests. The minimal requirement is that the code is stable for a 
range of simulated volumes and grid spacings a, such that an extrapolation to the continuum 
may be attempted (keeping the simulated volume fixed but sending a 0). Our code fulfills 
this property. 

Furthermore, one has to test whether this continuum extrapolation corresponds to the 
correct physical result in simple test cases. One such test case is provided by the 0-1-1 
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dimensional model discussed in section IIIBI Using initial conditions of uniform energy 
density in the 2+1 dimensional numerical code, the temperature evolution should match 
that of Eq. ([8]), for which it is straightforward to write an independent numerical solver. 
Our 2+1 dimensional code passes this test, for small and large ri/s and different values for 
Tn, Ai. 

The above test is non-trivial in the sense that it allows to check the implementation 
of nonlinearities in the hydrodynamic model. However, it does not probe the dynamics 
of the model, since e.g. all velocities are vanishing. Therefore, another test that one can 
(and should!) conduct is to study the dynamics of the model against that of linearized 



691 ] for similar considerations). 

but non- 



hydrodynamics (this test was first outlined in Ref. [49|; see 
More specifically, let us consider a viscous background "solution" with u 
vanishing e(r),n^(r) obeying Eq. ([8]). To first order in small fluctuations 5e,5u^,5Ii^^ 
around this background the set of equations become 



C%e + \drli\ + An| + (6 + p + Ui\)dr 

'c%e + \d^Il\ + + (e + p + \ll\)d^ 



6u^ + c%Se + dM" 



2t 

dr + 



6uy + cldy6e + di6Uy' 



1 + c 



3r Tu 



3r Tn 



6W - 



27] 

3rnr ' 4r] 



+ 



4?7 

3ttu 
6e 



6e + 



+ 



[e + p) + -n| 



T ^ 



4rn" 
2r] 



-n 



6e 



+ 



3r Tn 



3rn 3 « 

3rn 3 
V 



5Wy + ^ {djuy + dydu'') 









0, 



Tn 



(29) 



where we have put Ai = and assumed a constant for simplicity. Noting that 511^^ = 
SUl - (5n^^ from 6Ui/^ = 0, Eq. ([29D closed set of linear, but coupled differential 

equations for the fluctuations 6e, 6u^, Su^, 5n|, 611^^, 611^^. Doing a Fourier transform. 



6e{T,x,y) 



{2ny 



5e(r, r , F) 



(30) 



(and likewise for the other fluctuations), Eq. fl29|) constitute coupled ordinary differential 
equations for each mode doublet k = (A;^, A;^), which again are straightforward to solve with 
standard numerical methods [sH (and analytically for ideal hydrodynamics). 
A useful test observable is the correlation function 



/(r, Xi,X2) 



(5e(r, xi)(5e(r, xa)) 



e r 



(31) 



where () denotes an ensemble average over initial conditions 5e\^^^^. In particular, let us 
study initial conditions where 5e is given by Gaussian random noise with standard deviation 
A, 

/(ro,xi,X2) = AV(xi-X2) (32) 

and all other fluctuations vanish initially. These initial conditions are readily implemented 
both for the full 2+1 dimensional hydrodynamic code as well as for the linearized system 
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k [GeV] 



FIG. 3: (Color online) The correlation function /(r, k) as a function of momentum k = |k| for a 
lattice with a = 1 GeV~^, 64^ sites and averaged over 30 initial configurations (symbols), compared 
to the result from the linearized hydrodynamic equations (lines). 



Eq. (p9l) . As the system evolves to finite time r, both approaches have to give the same 
correlation function / as long as the linearized treatment is applicable, and hence Eq. (!29|) 
can be used to test the dynamics of the full numerical code. 

In practice, note that for the above construction / can only depend on the difference of 
coordinates, 



(5e(r, xi)5e(r, xa)) 



e r 



/(r,xi -xa) = j 



3«k-(xi-X2) 



/(^,k) 



and therefore in Fourier-space 

/(r,k)52(k') 



((5e(r,k)^e(r,k--k)) 
(27r)2e(r)2 



(33) 



(34) 



In the full 2-1-1 dimensional numerical code which is discretized on a space-time lattice, 
5^(k') is regular for any finite a, and one can maximize the signal for /(r, k) by calculating 
the r.h.s. of Eq. flM|) for k' = 0. Similarly, one solution 5e(r, k) per k mode is sufficient 
calculate /(t, k) for the linearized system Eq. fl29|) . 

The above initial conditions imply /(r = ro,k) = A^, but for finite times characteristic 
peaks develop as a function of |k|, whose position, height and width are sensitive to the values 
of c^, Tn, Tj/ s and of course the correct implementation of the hydrodynamic equations. The 
comparison between full numerics and linearized treatment shown in Fig. [3] suggests that 
our code also passes this test^. 

Finally, for the case of ideal hydrodynamics, analytic solutions to the hydrodynamic 
equations are known [t^, 71, 72|. Specifically, the code for central collisions has been 
found to agree with the results from Ref. [t^I for ideal hydrodynamics. Since our code agrees 
with Ref. in the case of central collisions and when dropping the appropriate terms in 
the equations ([6]), this provides yet another test on our numerics. 



^ Note that a small numerical error occurred in the linearized hydrodynamic solver and the corresponding 
figure in Ref. 22] . This error has been corrected in Fig. [31 
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To summarize, after conducting the above tests we are reasonably confident that our 
numerical 2+1 dimensional code solves the relativistic viscous hydrodynamic equations (fT4|) 
correctly. This completes the setup of a viscous hydrodynamic description of relativistic 
heavy-ion collisions. In the following sections, we will review comparisons of viscous hydro- 
dynamic simulations to experimental data, for both Glauber and CGC initial conditions. 



III. INITIAL CONDITIONS: GLAUBER MODEL VS. CGC 

A. The Glauber model 

In the Glauber model [7|, the starting point is the Woods-Saxon density distribution for 
nuclei, 

Pa{^) = -7-. 777^, p w p (35) 

1 +exp[(|x| -Ro)/x] 

where for a gold nucleus with weight A = 197 we use Rq = 6.4 fm and x = 0.54 fm. The 
parameter po is chosen such that / (i^xpyi(x) = A. One can then define the nuclear thickness 
function 

/oo 
dzpAi^) (36) 
-oo 

and subsequently the number density of nucleons participating in the collision (npart) and 
the number density of binary collisions (ncou)- For a collision of two nuclei with weight A 
at an impact parameter b, one has 



Ta 1^+2'^ 



gTa {x - 
1 



(jTa [x + 

' — 



ncon{x,y,b) 



aTA Ix + ^, yj Ta lx-^,y 



(37) 



where a is the nucleon-nucleon cross section. We assume a ~ 40 mb for Au+Au collisions 
at y/s = 200 GeV per nucleon pair. 

While the total number of participating nucleons A^part(&) = / dxdyn^artix^y.b) will be 
used to characterize the centrality class of the collision, as an initial condition for the energy 
density we will only use the parametrization 



e(r = To, X, y, b) = const x ncoii(a;, y, b), 



(38) 



since it gives a sensible description of the multiplicity distribution of experimental data, 
as will be discussed later on. In the following, "Glauber-model initial condition" is used 
synonymous to Eq. (1381) . 

The constant in Eq. (1381) is chosen such that the central energy density for zero impact 
parameter, e(r = tq, 0,0,0) corresponds to a predefined temperature Tj via the equation 
of state. This temperature will be treated as a free parameter and is eventually fixed by 
matching to experimental data on the multiplicity. 
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B. The CGC model 



The other model commonly used to obtain initial conditions for hydrodynamics is the so- 
called Color- Glass- Condensate approach, based on ideas of gluon saturation at high energies. 
In particular, we use a modified version of the KLN (Kharzeev-Levin-Nardi) /cr-factorization 
approach [tJ], due to Drescher et al. [t^]. We follow exactly the procedure described in [HI] 
and in fact we use the same numerical code, provided to us by the authors and only slightly 
modified to output initial conditions suitable for input into our viscous hydrodynamics 
program. 

In this model, the number density of gluons produced in a collision of two nuclei with 
atomic weight A is given by 

=^J^p d'kr asikr) M^u (Pt + ^t)V^; xt) M^2, (pt - kr^/A; x^) 

(39) 

where p-r and Y are the transverse momentum and rapidity of the produced gluons, respec- 
tively. Xi^2 = Pt ^ exp(±y)/y^ is the momentum fraction of the colliding gluon ladders 
with ^/s the center of mass collision energy and ^^(/cr) is the strong coupling constant at 
momentum scale = Ikrl- 

The value of the normalization constant Af is unimportant here, since as for Glauber 
initial conditions, we treat the overall normalization of the initial energy density distribution 
as a free parameter. The unintegrated gluon distribution functions are taken as 



as(Qs) niax((52, k? 



X, 4; x^) = -—— ^^-^ P(x^)(l - xr , (40) 



P(xt) is the probability of finding at least one nucleon at transverse position x^, taken from 
the definition for npart 



A 



PM = l-[l-^) , (41) 

where Ta and a are as defined in the previous section. 

The saturation scale at a given momentum fraction x and transverse coordinate x^ is 
given by 

The growth speed is taken to be A = 0.288. 

The initial conditions for hydrodynamic evolution require that we specify the energy 
density in the transverse plane at some initial proper time tq at which the medium has 
thermalized. Eq. (!39|) . on the other hand, is in principle valid at a time r<j = l/Qs at which 
the medium is likely not yet in thermal equilibrium. To obtain the desired initial conditions, 
we again follow [13] and assume that the number of gluons is effectively conserved during the 
evolution from Tg to tq and so the number density profile is the same at both times, scaled 
by the one-dimensional Bjorken expansion n(ro) = ^n(Ts). The energy density can then be 
obtained from the number density through thermodynamic relations — it is proportional to 
the number density to the 4/3 power. Again, we take the overall normalization as a free 
parameter, so the initial energy density is finally given as 



^{'^ = ^t, b) = const X 



dN, 



a 



[XT, h) 



4/3 



(43) 
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(a) (b) 




FIG. 4: (Color online) Left: The initial spatial anisotropy for the Glauber and CGC model. Right: 
The time evolution of the spatial and momentum anisotropy for a collision with & = 7 fm in ideal 
hydrodynamics. 



where the number density is given by Eq. (139|) evaluated at central rapidity Y = 0. 

As a final comment, it should be pointed out that the original version of the CGC, the 
McLerran-Venugopalan model 75|, l76|, differs from the KLN ansatz we used here, as will be 
discussed in the next-section. 



C. Spatial and momentum anisotropy 

One of the key parameters discussed in the following is the eccentricity (or spatial 
anisotropy) of the collision geometry. Following 0], we define it as 

e ^MH^ (44) 

where {)e denotes an averaging procedure over space with the energy density e as a weighting 
factor. Shown in Fig. H^, a plot of Cx for different centralities highlights the quantitative 
difference between the initial energy density from the Glauber and CGC model, Eq. (138|) 
and Eq. (H3!) . respectively. As can be seen from this figure, the CGC model generally gives 
a higher spatial anisotropy than the Glauber model. Note that the results for the CGC 
model shown here are extreme in the sense that the McLerran-Venugopalan model gives 
spatial eccentricities which essentially match the ones from the Glauber model This 
allows us to use the difference between the CGC and Glauber models as an indication of 
the systematic theoretical error stemming from the ignorance of the correct physical initial 
condition. 

Hydrodynamics converts pressure gradients into fluid velocities, and hence one expects 
the spatial anisotropy to decrease at the expense of a momentum anisotropy (which is related 
to the magnitude of the elliptic flow). We follow 7^] in defining a momentum anisotropy 
according to 
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FIG. 5: (Color online) Spatial and momentum anisotropy for the Glauber model at 6 = 7 fm with 
Ti = 0.353 GeV, tq = 1 fm/c and various values for the viscosity (grid spacing a = 2 GeV~^). 
(a): The dependence on the initialization value of the shear tensor: shown are results for vanishing 
initial value (n[^'^^ = 0) and Navier-Stokes initial value {Il'^^i^ 7^ 0), given in Eq. Q. (b): The 
dependence on the choice of value for rn, Ai: shown are results for rn = Ai = (labelled "IS") 
and rn = HiH^^2, Xi = ^ (labelled "AdS"). For rn = ^i2^^R, the results for Ai = (not 
shown) would be indistinguishable by bare eye from those for Ai = 



where we stress that here () denotes spatial averaging with weight factor unity. Fig.|3)D shows 
the time evolution in ideal hydrodynamics {rj/s -C 1) of both the spatial and momentum 
anisotropies for a heavy-ion collision at 6 = 7 fm modeled through Glauber and CGC 
initial conditions. As one can see, for the same impact parameter, the higher initial spatial 
anisotropy for the CGC model eventually leads to a higher momentum anisotropy than the 
Glauber model. Using a quasiparticle interpretation where the energy momentum tensor is 
given by 

the momentum anisotropy Cp can be approximately related to the integrated elliptic flow 



f2 (6), with a proportionality factor of ~ 2 [78|, l79|]. We find this proportionality to be 



maintained even for non- vanishing shear viscosity, as can be seen in Fig. [HI 



IV. RESULTS 

A. Which parameters matter? 

In the following, we will attempt to obtain limits on the mean value (throughout the 
hydrodynamic evolution) of the ratio rj/s from experimental data. While e.g. temperature 
variations of rj/s are to be expected in the real physical systems, probing for such variations 
would invariable force us to introduce more unknown parameters. We prefer to leave this 
program for future studies once robust results for the mean value of rj/s exist. Having fixed 
the equation of state and the freeze-out procedure as explained in the previous sections, the 
remaining choices that have to be made in the hydrodynamic model are the 
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Initial energy density profile: Glauber or CGC 

Initial value of shear tensor: vanishing or Navier-Stokes value 

Hydro dynamic starting time Tq 

Second-order coefficients: relaxation time rn and Ai 

Ansatz for non-equilibrium particle distribution Eq. ( fT9l) 



where it is to be understood that we fix the initial energy density normalization (Tj) and 
the freeze-out temperature Tf such that the model provides a reasonable description of the 
experimental data on multiplicity and < px >■ Historically, a strong emphasis has been 
put on requiring a small value of tq for ideal hydrodynamics [sil, 83] ■ For this reason, we 



will discuss the dependence on tq separately in section lIVDi A good indicator for which 
parameters matter is the momentum anisotropy since it is very sensitive to the value of rj/s. 
From Fig. [Hone therefore immediately concludes that the choice of Glauber or CGC initial 
conditions is important since it has a large effect on e^. Fortunately, most of the other 
choices turn out not to have a strong inffuence on the resulting V2 coefficient and hence the 
extracted rj/s. In the following we test for this sensitivity by studying for a "generic" 
heavy-ion collision of two gold nuclei, modeled by Glauber initial conditions at an initial 
starting temperature of Tj = 0.353, an impact parameter of 6 = 7 fm, and various choices 
of the above parameters. 

Fig. O shows the time evolution of e^jCp for various values of rj/s. From these plots, it 
can be seen that Cp (and hence V2) clearly is sensitive to the value of r]/s, suggesting that it 
can be a useful observable to determine the viscosity of the fluid from experiment. However, 
in order to be a useful probe of the fluid viscosity, the dependence of the flnal value of Cp 
on other parameters should be much weaker than the dependence on rj/s. In Fig. [5^ we 
show Cp, calculated for H^'^(ro) = and H^^(ro) equal to the Navier-Stokes value, Eq. (Q. 
As can be seen from this flgure, the resulting anisotropies are essentially independent of 
this choice, corroborating the flnding in Ref. jlo], 3l[. Similarly, in Fig. we show e„ 



calculated in simulations where the values of the second-order transport coefficients were 
either those of a weakly-coupled Miiller-Israel- Stewart theory (rn = 6^, Ai = 0) or those 
inspired by a strongly coupled M = 4 SYM plasma (rn = 2(2 — In 2)^, Ai = 2^)- Again, 
the dependence of Cp on the choice of the values of rn, Ai can be seen to be very weak for the 
values ofrj/s shown here. This result is in stark contrast to the findings of Ref. 2^, where a 



large sensitivity on the value of rn was found. However, recall that Ref. 29j] used evolution 
equations that differ from Eq. ([6]) and in particular do not respect conformal invariance. As 
argued in section IIIBl it is therefore expected to encounter anomalously large sensitivity to 
the value of the second order transport coefficients. 

To study the dependence of results on the ansatz of the non-equilibrium particle distri- 
bution function f[T^ . one would want to quantify the effect of neglecting terms of higher 
order in momenta in Eq. (flQl) . To estimate this, let us rewrite EcfN/cPp = EcPN^'^^ d^p + 
Ed^N^^'^ /d'^p, where N^^'^ contains only the equilibrium part where f{x^^p^) = fo (^^^), 
and perform a Pade-type resummation, 

^SjyPade ^ ^3^(0) ^ 
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FIG. 6: (Color online) Charged hadron elliptic flow for the Glauber model at h 
Ti = 0.353 GeV, tq = 1 fm/c and various viscosities. 



7 fm with 



Since Eq. ( 1471) contains powers of momenta to all orders when re-expanded, the difference 
between the ansatz (IT^ and the Pade resummed particle spectra can give a handle on the 
systematic error of the truncation used in Eq. (IT^ . Shown in Fig. El this difference suggests 
that this systematic error is small for momenta px ^ 2.5 GeV. Therefore, we do not expect 
our results to have a large systematic uncertainty coming from the particular ansatz f|T9|) 
for these momenta. 

To summarize, for values of rj/s ^ 0.2, the results for the momentum anisotropy are 
essentially insensitive to the choices for the second-order transport coefficients rn,Ai and 
the initialization of the shear tensor II'^'^[t = Tq). Conversely, is sensitive to the value 
of viscosity and the choice of initial energy density profile (initial eccentricity). Since the 
physical initial condition is currently unknown, this dependence will turn out to be the 
dominant systematic uncertainty in determining rj/s from experimental data. 



B. Multiplicity and radial flow 

As outlined in the introduction, we want to match the hydrodynamic model to experimen- 
tal data for the multiplicity, thereby fixing the constant in Eqs. f l55]) . fH51) . This translates 
to fixing an initial central temperature Tj for 6 = 0, which we will quote in the following. 

For a constant speed of sound, the evolution for ideal hydrodynamics is isentropic, while 
for viscous hydrodynamics additional entropy is produced. Since the multiplicity is a mea- 
sure of the entropy of the system, one expects an increase of multiplicity for viscous compared 
to ideal hydrodynamic evolution. This increase in final multiplicity has been measured as 
a function of rj/s for the semi-realistic speed of sound Fig. [T]in central heavy-ion collisions 
in Ref. [2H, and found to be approximately^ a factor of 0.751]/ s. (See Ref. [H, 83| for re- 
lated calculations in simplified models.) Reducing Tj accordingly therefore ensures that for 



^ The quoted fraction is for a hydrodynamic starting time of tq = 1 fm/c. Reducing tq leads to considerably 
larger entropy production. 
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FIG. 7: (Color online) Centrality dependence of total multiplicity dN/dY and < px > for 
7r+, 7r~, _fC+, K^,p and p from PHENIX [8J] for Au+Au collisions at ^ = 200 GeV, compared to 
the viscous hydrodynamic model and various rj/s, for Glauber initial conditions and CGC initial 
conditions. The model parameters used here are tq = 1 fm/c, rn = 6r//s, Ai = 0, = 140 MeV 
and adjusted Tj (see Table H]). 



viscous hydrodynamics, the multiplicity in central collisions will stay close to that of ideal 
hydrodynamics. 

Hydrodynamics gradually converts pressure gradients into flow velocities, which in turn 
relate to the mean particle momenta. Starting at a predefined time tq and requiring the 
hydrodynamic model spectra to match the experimental data on particle < pt > then fixes 
the freeze-out temperature Tf. 

For both Glauber-type and CGC-type model initial conditions, the experimental impact 
parameter dependence of the multiplicity and < px > is reasonably well parametrized for 
both ideal hydrodynamics as well as viscous hydrodynamics provided Tj is adjusted accord- 
ingly (see Fig. [7]). The values for Tj used in the simulations are compiled in Table[Tl We recall 
that no chemical potential is included in our equation of state, prohibiting a distinction be- 
tween particles and anti-particles, and chemical and kinetic freeze-out of particles occurs at 
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Initial condition 




Ti [GeV] 


Tf [GeV] 


To [fm/cj 


a [GeV ^] 


Glauber 




0.340 


0.14 


1 


2 


Glauber 


0.08 


0.333 


0.14 


1 


2 


Glauber 


0.16 


0.327 


0.14 


1 


2 


CGC 


10-4 


0.310 


0.14 


1 


2 


CGC 


0.08 


0.304 


0.14 


1 


2 


CGC 


0.16 


0.299 


0.14 


1 


2 


CGC 


0.24 


0.293 


0.14 


1 


2 



TABLE I: Summary of parameters used for the viscous hydrodynamics simulations 

the same temperature. Furthermore, approximating the equihbrium particle-distributions 
for bosons by a Boltzmann distribution (|T9l) leads to small, but consistent underestimation 
of the multiplicity of light particles, such as pions. For these reasons, it does not make sense 
to attempt a precision fit to experimental data, especially for pions and protons. Rather, 
we have aimed for a sensible description of the overall centrality dependence of multiplicity 
and < pt > of kaons. 

Note that in particular for the CGC model one could achieve a better fit to the data on 
mean < > by increasing the freeze-out temperature by ~ 10 MeV. This would also lead 
to a decrease in elliptic flow for this model. However, to facilitate comparison between the 
CGC and Glauber initial conditions, we have kept Tf the same for both models. 



C. Elliptic flow 

Having fixed the parameters TQ,Ti, Tf for a given rj/s to provide a reasonable description 
of the experimental data, a sensible comparison between the model and experimental results 
for the elliptic fiow coefficient can be attempted. For charged hadrons, the integrated and 
minimum-bias f 2 coefficients are shown in Fig. [H] for Glauber and CGC initial conditions. 
As noted in section IIII C[ charged hadron f 2°* turns out to be very well reproduced by the 
momentum eccentricity | Cp, evaluated when the last fiuid cell has cooled below Tf. This 
agreement is independent from impact parameter or viscosity and hence may serve as a more 
direct method on obtaining an estimate for f ™* if one cannot (or does not want to) make 
use of the Cooper- Frye freeze-out procedure described in section III G[ 

The comparison of the hydrodynamic model to experimental data with 90% confidence 
level systematic error bars from PHOBOS jssi for the integrated elliptic fiow in Fig. [S] 
suggests a maximum value of 77/5 ~ 0.16 for Glauber-type and 77/5 ~ 0.24 for CGC-type 
initial conditions. Whereas for Glauber initial conditions, ideal hydrodynamics {r]/s ~ 0) 
gives results consistent with PHOBOS data, for CGC initial conditions zero viscosity does 
not give a good fit to the data, which is consistent with previous findings ^]. 

For minimum-bias V2, to date only experimental data using the event-plane method are 
available, where the statistical, but not the systematic error of that measurement is directly 
accessible. The dominant source of systematic error is associated with the presence of so- 



called non-fiow effects [86|. Recent results from STAR suggest that removal of these non-flow 



effects imply a reduction of the event-plane minimum bias V2 by 20 percent [87|, l88|. For 



charged hadrons, a comparison of both the event-plane and the estimated non-flow corrected 
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FIG. 8: (Color online) Comparison of hydrodynamic models to experimental data on charged 
hadron integrated (left) and minimum bias (right) elliptic flow by PHOBOS [85|] and STAR [83], 
respectively. STAR event plane data has been reduced by 20 percent to estimate the removal 



of non-flow contributions 



87|, 



The line thickness for the hydrodynamic model curves is an 
estimate of the accumulated numerical error (due to, e.g., finite grid spacing). The integrated V2 
coefficient from the hydrodynamic models (full lines) is well reproduced by ^e^ (dots); indeed, the 
difference between the full lines and dots gives an estimate of the systematic uncertainty of the 
freeze-out prescription. 



experimental data from STAR with the hydrodynamic model is shown in Fig. [8l 

For Glauber-type initial conditions, the data on minimum-bias V2 for charged hadrons 
is consistent with the hydrodynamic model for viscosities in the range rijs G [0,0.1], while 
for the CGC case the respective range is t^/s G [0.08,0.2]. It is interesting to note that 
for Glauber-type initial conditions, experimental data for both the integrated as well as the 
minimum-bias elliptic flow coefficient (corrected for non-flow effects) seem to be reproduced 
best^ by a hydrodynamic model with rj/ s = 0.08 ~ This number has first appeared in the 



^ In Ref. a lower value of 77/s for the Glauber model was reported. The results for viscous hydrodynamics 
shown in Fig. [5] are identical to Ref. 2^, but the new STAR data with non-flow corrections became 
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FIG. 9: (Color online) Momentum anisotropy (a) and elliptic flow for charged hadrons (b) for 
6 = 7 fm, r}/ s = 0.08 and different hydrodynamic initialization times tq. Horizontal light gray lines 
in (a) are visual aids to compare the final value of Cp. As can be seen from these plots, neither 
the final Cp nor the charged hadron V2 depend sensitively on the value of tq if the same energy 
distribution is used as initial condition at the respective initialization times. Simulation parameters 
were % = 0.29 GeV, Tj = 0.14 GeV for tq = 2 fm/c, Ti = 0.36 GeV, T/ = 0.15 GeV for tq = 1 
fm/c, and Ti = 0.43 GeV, T/ = 0.16 GeV for tq = 0.5 fm/c. 



gauge/string duality context [10| and has been conjectured to be the universal lower bound 
on r]/s for any quantum field theory at finite temperature and zero chemical potential (soj ]. 
For CGC-type initial conditions, the charged hadron V2 data seems to favor a hydrodynamic 
model with t]/s ^ 0.16, well above this bound. 



D. Early vs. late thermalization 

Currently, there seems to be a common misunderstanding in the heavy-ion community 
that hydrodynamic models can universally only reproduce experimental data if they are ini- 
tialized at early times tq < 1 fm/c. This notion has been labeled "early thermalization" and 
continues to create a lot of confusion. In this section, we argue that the matching of hydro- 
dynamics to data itself does not require Tq < 1 fm/c. It is the additional assumptions about 
pre-equilibrium dynamics that lead to this conclusion for the Glauber initial conditions. 

Performing hydrodynamic simulations in the way we have described earlier, the energy 
density distribution is specified by either the Glauber or CGC model at an initial time tq. 
In Fig. M we show the result for the elliptic flow coefficient (or the momentum anisotropy) 
for three different values of tq, namely 0.5, 1 and 2 fm/c, where also Tj and Tf have been 
changed in order to obtain roughly the same multiplicity and mean pt for each tq. As can 
be seen from this figure, the resulting final elliptic fiow coefficient is essentially independent 
of the choice of tq. In particular, this implies that experimental data for bulk quantities 
can be reproduced by hydrodynamic models also for large initialization times, so no early 
thermalization assumption is needed. 



available only after 



22l | had been published. 
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However, it is true that the above procedure assumes that the energy density distribu- 
tion remains unchanged up to the starting time of hydrodynamics, which arguably becomes 
increasingly inaccurate for larger tq. It has therefore been suggested |81|] to mimic the 
pre-hydro time evolution of the energy density distribution by assuming free-streaming of 
partons. Assuming free-streaming gives the maximal contrast to assuming hydrodynamic 
evolution, since the latter corresponds to very strong interactions while the former corre- 
sponds to no parton interactions at all. Indeed, one can calculate the effect of the free- 



streaming evolution on the spatial anisotropy, finding [81 



This implies that the spatial anisotropy decreases with time, whereas one can show that 
free-streaming does not lead to a build-up of Cp. In other words, the eccentricity gets diluted 
without producing elliptic flow, such that once hydrodynamic evolution starts, it will not 
lead to as much V2 as it would have without the dilution effect*. It is tempting to conclude 
from this that by comparing to experimental data on elliptic flow one could place an upper 
bound on the maximally allowed dilution time, and interpret this as the thermalization time 
of the system. One should be aware, however, that this bound will depend on the assumption 
made about the pre-hydro evolution. Furthermore, one should take into account the fact 
that the initial state of the system remains unknown. For instance, the system could start 
with an energy density distribution similar to the CGC model, which has a fairly large 
eccentricity. Fig. [10] shows that when allowing the eccentricity to get diluted according to 
Eq. (jUD, it takes a time of r ~ 1.5 fm/c until the eccentricity has shrunk to that of the 
Glauber model. This implies that even when assuming no particle interactions (no elliptic 
flow build-up) for the first stage of the system evolution, one can get eccentricities which are 
Glauber-like after waiting for a significant fraction of the system life time. Allowing at least 
some particle interactions (which is probably more realistic), one expects some build-up of 
elliptic flow already in the dilution (or pre-equilibrium) phase, and therefore dilution (or 
"thermalization") times of r ~ 2 fm/c seem not to be incompatible with the observed final 
elliptic flow even for non- vanishing viscosity. 



V. SUMMARY AND CONCLUSIONS 

In this article, we applied conformal relativistic viscous hydrodynamics to simulate 
Au+Au collisions at RHIC at energies of \fs = 200 GeV per nucleon pair. Besides one 
first-order transport coefficient (the shear viscosity) in general there are five second-order 
transport coefficients in this theory, for which one would have to supply values. We provided 
arguments that physical observables in the parameter range accessible to hydrodynamics 
(low momenta, central to semi-central collisions) do not seem to be strongly dependent on 



It seems that if one forces the energy-momentum tensor at the end of free-streaming period to match to 
that of ideal hydrodynamics (instantaneous thermahzatiqn), the resuhing fluid velocities are anisotropic, 



i.e. correspond to a non- vanishing elliptic flow coefficient [89|,|90|. It is possible that this effect stems from 
neglecting velocity gradients (viscous hydrodynamic corrections) in the matching process. We ignore the 
complications of the detailed matching from free-streaming to hydrodynamics in the following. 
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FIG. 10: (Color online) Spatial eccentricity for the Glauber and CGC model compared to evolving 
the CGC model according to Eq. (|48p for r = 1.5 fm/c. This implies that starting with Glauber-type 
initial conditions at tq > 1 fm/c may not be unreasonable. 

specific (reasonable) choices for any of these second-order coefficients. On the other hand, 
we do find a pronounced dependence of the eUiptic flow coefficient on the ratio of shear vis- 
cosity over entropy density, which suggests that by combining viscous hydrodynamics and 
experimental data a measurement of the quark-gluon plasma viscosity may not be futile. 
However, we have shown that our ignorance about the precise distribution of energy density 
at the earliest stages of a heavy-ion collision introduces a large systematic uncertainty in 
the final elliptic flow of the hydrodynamic model. Adding to this is the considerable ex- 
perimental uncertainty pertaining to the removing of non-flow contributions to the elliptic 
flow. For these reasons, we are unable to make precise statements about the value of the 
shear viscosity of the quark-gluon plasma and in particular cannot place a firm lower bound 
on Tj/s. Indeed, our hydrodynamic models seem to be able to consistently describe experi- 
mental data for multiplicity, radial flow and elliptic flow of bulk charged hadrons for a wide 
range of viscosity over entropy ratios, 

- = 0.1 ± 0.1 (theory) ± 0.08 (experiment), (49) 

where we estimated the systematic uncertainties for both theory and experiment from the 
results shown in Fig. [Hi We stress that Eq. fH^ does not account for physics not included in 
our model, such as finite chemical potential, bulk viscosity, heat flow, hadron cascades, three- 
dimensional fluid dynamic effects and possibly many more. Consistent inclusion of all these 
may result in changes of the central value and theory uncertainty in Eq. (H9|) . Nevertheless, 
none of the mentioned refinements is currently expected to dramatically increase the elliptic 
flow coefficient (though some increase may be expected when e.g. implementing partial 
chemical equilibrium [9l|). Therefore, we seem to be able to exclude viscosities of r^/s ^ 0.5 
with high confidence, which indicates that the quark-gluon plasma displays less friction than 



931,194,195 



any other known laboratory fluid [80|, l92j. Other groups have come to similar conclusions 



To better quantify the shear viscosity of the quark-gluon plasma at RHIC calls for more 
work, both in theory and experiment. On the theory side, a promising route seems to be the 
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stud y of fluctuations and comparing to existing experimental data |85l.l93. [gel. [97l.l98l.l99l. [Tool. 



Toil . For instance, it might be interesting to investigate the critical value oirj/ s for the onset 



of turbulen ce in heavy-ion collisions and explore possible consequences of fully developed 



turbulence [102| . However, maybe most importantly, a more thorough understanding of 
the earliest stages of a heavy-ion collision, in particular thermalization, could fix the initial 
conditions for hydrodynamics and hence dramatically reduce the theoretical uncertainty in 
final observables. 

Leaving these ideas for future work, we stress that with the advent of conformal relativis- 
tic viscous hydrodynamics at least the uncertainties of the hydrodynamic evolution itself 
now seem to be under control. We hope that this serves as another step towards a better 
understanding of the dynamics of relativistic heavy-ion collisions. 



Note on version 4 of preprint 

An error was found in the freeze out routine used to obtain the original results (see the 
pubhshed erratum [Phys. Rev. C 79, 039903(E) (2009)]). The error was corrected and 
the matching to the experimental data was redone. The only changes are the values of the 
parameters listed in Table [U and slightly revised Figs. El [71 [HI and [9b, which have been 
updated in this version 4. 
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